Eutectic colony formation: A phase field study 
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Eutectic two-phase cells, also known as eutectic colonies, are commonly observed during the so- 
lidification of ternary alloys when the composition is close to a binary eutectic valley. In analogy 
with the solidification cells formed in dilute binary alloys, colony formation is triggered by a mor- 
phological instability of a macroscopically planar eutectic solidification front due to the rejection 
by both solid phases of a ternary impurity that diffuses in the liquid. Here we develop a phase- 
field model of a binary eutectic with a dilute ternary impurity and we investigate by dynamical 
simulations both the initial linear regime of this instability, and the subsequent highly nonlinear 
evolution of the interface that leads to fully developed two-phase cells with a spacing much larger 
than the lamellar spacing. We find a good overall agreement with our recent linear stability analysis 
[M. Plapp and A. Karma, Phys. Rev. E 60, 6865 (1999)], which predicts a destabilization of the 
front by long-wavelength modes that may be stationary or oscillatory. A fine comparison, however, 
reveals that the assumption commonly attributed to Cahn that lamella grow perpendicular to the 
envelope of the solidification front is weakly violated in the phase-field simulations. We show that, 
even though weak, this violation has an important quantitative effect on the stability properties of 
the eutectic front. We also investigate the dynamics of fully developed colonies and find that the 
large-scale envelope of the composite eutectic front does not converge to a steady state, but exhibits 
cell elimination and tip-splitting events up to the largest times simulated. 



I. INTRODUCTION 

Eutectic alloys can form a wealth of different two-phase 
patterns during solidification. These alloys are of inter- 
est to metallurgists jj| because of their low melting points 
and of the superior mechanical properties associated with 
a fine-scale composite microstructure. Moreover, eutec- 
tic growth is a non-trivial example of pattern formation 
outside of equilibrium that has attracted the attention of 
physicists over the last two decades. 

When the two solid phases (a and 0) of a binary eu- 
tectic alloy have rough interfaces with the liquid, solidifi- 
cation at or near the eutectic composition typically pro- 
duces a spatially periodic array structure consisting of 
lamellar plates of the two phases, or of rods of the phase 
with the smaller volume fraction embedded inside the 
matrix of the other phase. Since the pioneering mathe- 
matical analyses by Hillert and Jackson and Hunt [3|, 
which built on earlier work by Brandt Q and Zener j 5 j , 
it is well established that these lamellar and rod mor- 
phologies can grow cooperatively in steady state for a 
continuous range of eutectic spacings, with both phases 
helping each other to grow via the diffusive transport 
of the two chemical components in the liquid (coupled 
growth) . 

In directional solidification experiments, a sample con- 
taining the alloy is pulled at a constant velocity v p in 
an externally imposed temperature gradient of magni- 
tude G. In such experiments, coupled growth typically 
produces a composite front that is perpendicular to the 




FIG. 1. Eutectic colonies in a thin sample of the transpar- 
ent organic eutectic alloy CBr4 — C2CI6, doped with a small 
amount of the ternary impurity naphtalene (from Ref. jlil). 



temperature gradient, and planar on a scale much larger 
than the lamellar spacing A (defined as the width of 
the basic spatially repeating unit consisting of one a- 
lamella and one /3-lamella) . Analytical Q and numerical 
studies of the morphological stability of lamellar eu- 
tectics, as well as detailed experiments in a transparent 
organic system [§), have shown that the stable range of 
lamellar spacings is restricted by a long-wavelength in- 
stability leading to local lamellar termination at small 
A, and short-wavelength oscillatory instabilities at large 
A. These studies clearly demonstrate that a large-scale 
morphological instability of the composite front does not 
occur in a binary eutectic alloy. 

This picture is consistent with the experimental obser- 
vation that such a morphological instability occurs only 
when a small quantity of a ternary impurity is present, 
and when v p exceeds a critical value [p|-|l9|]. In a nonlin- 
ear regime, this instability gives rise to the formation of 
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two-phase solidification cells, also called eutectic colonies, 
with a typical width much larger than A. A typical ex- 
ample of such cells is shown in Fig. [|. 

Experimental measurements to date have con- 

sistently shown that the onset of colony formation can 
be relatively well described by a simple constitutional 
supercooling criterion with respect to the ternary impu- 
rity which predicts that instability occurs when 
G/vp falls below a critical value. This suggests that this 
instability may be qualitatively similar to the well-known 
Mullins-Sekerka instability of a monophase front during 
directional solidification of a dilute binary alloy [ p2[ . In a 
recent linear stability analysis of a sharp interface model 
p3| , however, we showed that the morphological insta- 
bility of a composite front in the presence of a dilute 
ternary impurity differs fundamentally from the instabil- 
ity of a monophase front, even though the onset of both 
instabilities is well predicted by constitutional supercool- 
ing. This analysis was based on the same procedure used 
previously by Datye and Langer || to analyze the stabil- 
ity of binary lamellar eutectics, where the basic degrees 
of freedom are the coordinates of the a-/3-liquid trijunc- 
tions. Our main finding was that the amplification of 
linear perturbations of the composite front can be either 
steady or oscillatory for experimentally relevant control 
parameters, in contrast to the classical Mullins-Sekerka 
instability where finite- wavelength perturbations are am- 
plified in a non-oscillatory way. 

Furthermore, in Ref. p3| ], we developed an "effective 
monophase front" formulation of the dynamics of the 
composite interface that shed light on the origin of this 
difference. We showed that the long-wavelength dynam- 
ics of the envelope of the composite front is governed by a 
free-boundary problem with boundary conditions for the 
concentration of the diffusing ternary impurity on the ef- 
fective front that can be obtained by averaging over the 
properties of the two solid phases. As a self-consistency 
check, we also showed that, when the wavelength of the 
perturbation is much larger than A, the linear stability 
analysis of this free-boundary problem gives identical re- 
sults to the full stability calculation expressed in terms 
of the trijunction coordinates. 

Not surprisingly, this free-boundary problem turns 
out to be very similar to the one governing a "true" 
monophase front in a dilute binary alloy. The non-trivial 
difference, however, is that the local lamellar spacing, 
which appears in the boundary condition for the ternary 
impurity on the front, constitutes an additional "internal 
degree of freedom" of the front that modifies its stability 
properties, and gives rise to the oscillatory modes. Phys- 
ically, this reflects the fact that the local temperature 
of the front depends on the local lamellar spacing A and 
that, in turn, the time rate of change of A depends on the 
shape of the front because of the geometrical constraints 
imposed by the equilibrium conditions for the angles be- 
tween interfaces at the trijunctions (Young's conditions). 

In a recent experimental study of a transparent organic 
model alloy, oscillatory patterns compatible with the re- 



sults of our linear stability analysis were indeed observed 
|Q. The same study, however, also revealed a wealth of 
other possible structures that can be associated with the 
instability of a planar front, and in particular localized 
two-phase fingers that may appear in an early stage of 
the morphological instability. 

The two main goals of the present study are to check 
the validity of our previous linear stability analysis |2lf | 
by direct simulation of the fundamental equations of mo- 
tion, and to investigate the nonlinear regime of colony 
formation. For this purpose, we develop a phase-field 
model for the directional solidification of a eutectic alloy 
with a dilute ternary impurity. Simulations of this model 
enable us to characterize quantitatively the amplification 
and decay of linear perturbations of the composite front 
and to study the complex interface dynamics leading to 
the formation of well-developed colonies. 

The phase-field method is by now a well-established 
technique for simulating solidification patterns [^i] p9| . 
In particular, it has already been applied to the investi- 
gation of multiphase solidification in eutectic and peri- 
tectic alloys pp|-p5| . The advantage of this method with 
respect to the boundary integral formalism used previ- 
ously to perform detailed simulations of eutectic growth 
structures is that ternary impurities can easily be in- 
cluded. Furthermore, the phase field method is able to 
handle automatically dramatic changes in the interface 
morphology such as lamella termination and creation, 
which are difficult to implement in the boundary inte- 
gral approach. 

The phase-field model presented in this paper is specif- 
ically designed for computational efficiency and therefore 
makes some simplifying assumptions. In particular, we 
use a generic eutectic phase diagram that is symmetric 
with respect to the exchange of the two solid phases, and 
we neglect crystallographic effects such as the anisotropy 
of the solid-liquid and solid-solid interfacial energies. The 
computational effort required to simulate fully developed 
colonies is nonetheless considerable since the two-phase 
cell spacing is one order of magnitude larger than A. For 
this reason, the largest simulations of such structures 
were carried out on a multi-processor CRAY T3E and 
took the equivalent of a few thousand hours of single- 
processor workstation time. 

The simulation results are found to be in good overall 
agreement with our sharp-interface linear stability anal- 
ysis for compositions close to the eutectic point, where 
the two solid phases have approximately equal volume 
fractions. We observe, indeed, the predicted large-scale 
oscillatory structures. Quantitatively, however, the sim- 
ulated growth rates differ from the predicted ones. A 
careful analysis of our simulation results, extrapolated 
to the limit of vanishing thickness of the diffuse inter- 
faces, allows us to pinpoint the origin of this discrepancy. 
In particular, our stability analysis uses the assumption 
that the lamellae grow normal to the large-scale growth 
front. This assumption is commonly attributed to Calm 
and was also used previously by Datye and Langer || for 



2 



their linear stability analysis of lamellar eutectics in bi- 
nary alloys. We find that, in our simulations, this rule is 
slightly violated. Hence, the stability analysis correctly 
describes all the qualitative features of the instability, but 
would have to be extended to include this effect in order 
to become quantitatively accurate. This violation also 
has important consequences for the stability of binary 
eutectics that will be discussed elsewhere. 

The linear instability of the planar front is followed by 
a nonlinear transient that leads to the formation of fully 
developed colonies. The nature of the transient depends 
on the composition. In simulations carried out at the 
eutectic point, the long- wavelength modes grow until the 
front becomes wavy and the first lamella terminations 
occur in the concave parts. Subsequently, the grooves 
deepen and the tips grow ahead of the front, such that 
the initial wavelength of the colonies corresponds to the 
linear mode that dominates the stability spectrum. In 
contrast, for off-eutectic compositions, the linear regime 
is much shorter, and localized two-phase fingers centered 
around a thin lamella of the minority phase grow rapidly 
ahead of the front and develop into colonies later on. 

Finally, once formed, the colonies have a quite well- 
defined average size and shape at both eutectic and off- 
eutectic compositions. However, the front does not settle 
down into a true steady state, but exhibits tip-splitting 
and cell elimination events, not unlike the monophase 
front of a dilute alloy in the absence of interfacial 
anisotropy [ ^6|j37[ . 

The remainder of this paper is organized as follows. 
In the next section, we introduce the phase-field model 
and analyze its sharp-interface limit. In section III, we 
present simulation results for stable steady-state lamel- 
lar growth that are used to test our model. Section IV 
contains a brief review of our sharp-interface linear sta- 
bility analysis p3|] , and a detailed comparison between 
the analytical and numerical results concerning the lin- 
ear stability of the eutectic front. Section V is devoted 
to the simulations of well-developed colony structures in 
a nonlinear regime. Finally, conclusions and an outlook 
for future work are given in section VI. 



II. PHASE-FIELD EQUATIONS AND 
SHARP-INTERFACE LIMIT 

We consider directional solidification of thin samples, 
as used in many experimental studies of pattern forma- 
tion during solidification |p|, ^2| , p^[ .This allows us to treat 
the problem as essentially two-dimensional and to ne- 
glect convection in the liquid. Furthermore, we assume 
that the rejection of latent heat during solidification does 
not appreciably modify the temperature field created by 
the experimental setup (frozen temperature approxima- 
tion), and hence that growth is limited by diffusion of 
the chemical constituents. 

We are interested in the behaviour of a ternary al- 
loy close to a binary eutectic trough in the phase dia- 



gram. Specifically, we will consider a very low concentra- 
tion of the third component, which can then be regarded 
as a dilute impurity. This allows us to neglect various 
cross-coupling terms between the ternary impurity and 
the components of the binary eutectic. In addition, we 
are more interested in generic aspects of two-phase cell 
formation than in modeling a specific material. Hence, 
we study a model eutectic alloy that has a symmetric 
phase diagram. This simplifies the setup of the phase- 
field model. 

The principles of the phase-field method have been de- 
scribed in detail in numerous publications p^-p5[. The 
idea is to distinguish between the different thermody- 
namic phases with the help of one or several scalar fields, 
the phase fields, that have fixed values in the bulk phases 
and vary continuously across smooth and diffuse inter- 
faces. A free energy functional suitable for the problem 
at hand is then constructed, and the equations of mo- 
tion for the fields are written in variational form. By 
now, various phase- field models for alloy solidification are 
available ||^,|2^-^|. In particular, much effort was spent 
to develop a thermodynamically consistent approach and 
to base the free energy functional on ideal or regular so- 
lution models [^7 32 34]]. In contrast, we are interested 
here mainly in the phase-field model as a computational 
tool. We will therefore use a strongly simplified model 
that is chosen for its computational efficiency, with the 
minimum of ingredients necessary to reproduce the main 
features of eutectic solidification with a ternary impu- 
rity. The parameters of the model are related to physical 
quantities by performing a sharp-interface limit. 

We choose as the set of dynamical field variables the 
concentration (in molar fraction) c(x, z, t) of one of the 
components of the binary eutectic, the concentration 
c{x,z,t) of impurities, and a single phase field 4>{x,z,t) 
that distinguishes between solid and liquid. To simplify 
the construction of the free energy functional, we define 
the scaled concentration u by 



u(x, z, t) 



c(x,z,t) - c E 
(c/3 - c a ) /2 



(1) 



where ce, c Q , and cp are the compositions of the liquid 
and the two solid phases in the pure binary eutectic at 
the eutectic temperature Te fl38[ . For a symmetric phase 
diagram, the scaled compositions of the two solids at Te 
are u = ±1. 

Building on a previous phase-field model for a binary 
eutectic fl30|j, we take the (dimensionless) free energy 
functional [t59l of the form 



F= dV 
'v 



W 2 o Wl , 

_n {Vu ) 2 + ^(V0) 2 



f((j>,u,c,T) 



(2) 



where V is the volume of the two-phase system. The di- 
mensionless free energy density /(</>, u, c, T) must have 



3 



three local minima to account for the three possible 
phases (liquid, a solid, and (3 solid), separated by po- 
tential barriers. We use the phase field to distinguish 
between solid and liquid, and the scaled concentration 
field to distinguish between the two solids. The gradient 
terms force the fields to vary continuously between the 
bulk equilibrium values and hence create interfaces of a 
characteristic thickness of order W u (solid-solid interface) 
and W<f, (solid- liquid interfaces) . In general, there should 
also be a gradient term for the ternary impurity. How- 
ever, we may omit this term for simplicity since c has 
no indicator function, but is slaved to the other fields; 
that is, for specified phase field (f>, concentration u and 
temperature T, the equilibrium value of c is known. 
A convenient choice for the free energy density is 



f[<p,u,c,T) = -— + — + f soX {u,c,T) 



i - h{4>) 



fu q (u,c,T). 



(3) 



Here, / so i and /n q are the bulk free energy densities in 
the solid and the liquid, respectively, and 



K4>) 



(4) 



is an interpolation function. The first two terms in Eq. 
(|J) generate a double well potential for (j) with minima 
at cj) = ±1. Since h(±l) = ±1, /(l, u, c, T) = f sol (u, c, T) 
and /(— 1, u, c, T) — fn q (u, c, T), such that (f> = +1 corre- 
sponds to the solid and <f> = — 1 to the liquid. Moreover, 
since /i'(±l) = 0, the equilibrium values of <f>, given by 
the solutions of df /d<j> = 0, always remain at 4> = ±1, 
independently of the values of / so i and fu q . 
For / so i and /n q , we take 



fn q (u,c, T) =u 2 /2 + bclnc- e t c, 
f so \(u,c, T) — a (u 2 - l) 2 + 6clnc 



e s c ■ 



(5) 

aAT/T E , (6) 



where AT = Te — T is the undercooling with respect to 
the binary eutectic point, and a, b, e s , e;, and a are con- 
stants that will be related to physical parameters by the 
construction of the phase diagram. This choice is moti- 
vated by the following considerations. Since there are two 
solid phases, / so i must have a double-well structure in u; 
in contrast, fn q has a single well. At the eutectic tem- 
perature and without impurities (c t o0), all three phases 
must have the same free energy; for T > Te (T < Te), 
the liquid minimum must be below (above) the solid min- 
ima. This is conveniently achieved by the last term on 
the right hand side of Eq. (Q) that simply shifts / so i with 
respect to /u q ; formally, a is equivalent to the latent heat. 

The impurity terms have a form that is equivalent to 
the dilute limit of a regular solution model. Indeed, the 
terms containing c correspond to the dilute approxima- 
tions for the entropy of mixing and the energy cost of the 
impurities, respectively, with e v representing the differ- 
ence in bond energies upon replacing a "solvent" atom 



by an impurity in phase v. The constant b, which sets 
the energy scale, should formally be proportional to the 
temperature. Since we are, however, only interested in a 
narrow temperature range around Tg, we simply use a 
constant. 

The various coefficients can be related to physical 
quantities through the construction of the phase diagram. 
The conditions for thermodynamical equilibrium between 
two distinct phases are (i) equal chemical potentials for 
the eutectic components, (ii) equal chemical potentials 
for the ternary impurity, and (iii) equal grand potential, 
i.e. 

fi s = df so i/du\ Us = m = dfu q /du\ ai , (7) 
A s = df so i/dc\~ s =iii= <9/n q /(9c| £! , (8) 

= /sol - HsU s - L~l s C s =0,1 = /liq - [llUl - /i/Q, (9) 

where u v and c u , v = 8,1, denote the equilibrium con- 
centrations in solid and liquid. These conditions can be 
geometrically described as a "common tangent plane" 
to the free energy surface, analogous to the well-known 
double-tangent construction for binary alloys. From Eq. 
(||), we get at once the standard partition relation for a 
dilute alloy, 



c s = Kcu 
with a partition coefficient given by 
K = exp[-(e s - e t )/b]. 



(10) 



(11) 



Next, from Eq. (R), under the assumption AT/Tg <C 1 
(i.e. for temperatures close to the eutectic temperature), 
we have iijCl and u s ~ ±1, and we get 



u s = — ± 1, 
8a 



(12) 



where the two signs correspond to the two distinct solid- 
liquid equilibria. Finally, using Eqs. @, (|l0|), and JT^), 
we obtain 

aAT/T E - 6(1 - K)ci ± u,. (13) 

Using the definition of AT, this can be rewritten as 

T = T E ±m u ui - rhci, (14) 

where m u and rh are the magnitudes of the liquidus 
slopes for the eutectic components and the impurity, re- 
spectively, 



dT 



dui 
dT 



dc. 



Te 
a 

6(1 - K)T E 
a 



(15) 
(16) 



Note that the scaled liquidus slope m u can be related to 
the "true" liquidus slope m in the phase diagram with 
the help of Eq. (g), 
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m(ca - c a )/2. 



(17) 



The parameter a controls the ratio of the liquidus and 
solidus slopes in the eutectic phase diagram; for simplic- 
ity, we will fix in the following a = 1/8, which gives a 
concentration jump across the interface that is indepen- 
dent of temperature (parallel liquidus and solidus lines). 
The parameter b, together with the partition coefficient 
K, fixes the ratio of eutectic and impurity liquidus slopes, 
m/m u = 6(1 — K). 

The equations of motion for the three fields are 



rd t (p(x,z,t) 



SF 



[x,z,ty 



d t u(x, z, t) = V M{<j>, u, c)V 



SF 



d t c(x,z,t) = V M(<j),u,c)V — 



Su(x, z,t) J 
SF 



Sc(x, z, t) 



(18) 
(19) 
(20) 



where SF/S- denotes the functional derivative with re- 
spect to the field •, r is a (microscopic) relaxation time, 
and M and M are the mobilities of the eutectic compo- 
nent and the ternary impurity, respectively. These vari- 
ational forms reflect the fact that the two concentrations 
are conserved fields, whereas the phase field can be seen 
as a non-conserved order parameter. The non-conserved 
phase-field simply relaxes towards its local equilibrium 
value. Indeed, by inserting Eq. (||) into Eq. ( |is| ) we 
obtain 

Td t 4> = W^V 2 + <j>/2 - 4> 3 /2 + h'(cf>)(h q - /soi). (21) 

The last term on the right hand side always drives the 
phase field to the value that corresponds to the lower 
local free energy density (we recall that h' > and that 
0=1 corresponds to the solid). 

The definition of the model is completed by the specifi- 
cation of the mobility functions M(<j>, u, c) and M(</>, u, c). 
The dependence of M and M on the phase field and the 
compositions allows us to obtain the desired diffusivities 
in the bulk phases. We want to simulate a one-sided 
model (i.e. vanishing diffusivity in the solid) with con- 
stant diffusivities for eutectic components and impurities 
in the liquid. This can be achieved by choosing 



M(4>,u,c) = D 
M(cj),u 7 c) = D 



1 - 



1 - 



1 + 



(22a) 
(22b) 



where D and D are the diffusion constants. Indeed, from 
the equations of motion, we get that in the liquid (0 = 

-1) 



AI 



d 2 fu q 
du 2 



Vu + VK 2 V(V 2 u) 



d t c = V [M 



(23) 
(24) 



In the first equation, we can neglect the term W 2 V(V 2 u) 
in the brackets on the right-hand side, since the diffusion 
pattern forms on a scale much larger than W u , and hence 
this term is small compared to (d 2 fu q /dc 2 )Wu. Using 
the expressions for the mobilities and /i; q , we obtain the 
desired result in the liquid, 



d t u = DV 2 u, 
dtc= DV 2 c. 



(25) 
(26) 



The exponent n in the mobility plays a role only in the 
interfacial region where the phase field varies, and chang- 
ing its value modifies the interface kinetics. This will be 
discussed in more detail below. 

When the thickness of the diffuse interfaces is much 
smaller than all other physical length scales, and in par- 
ticular the lamellar spacing A, the above phase-field equa- 
tions can be related to the more conventional sharp- 
interface equations of the macroscopic models of solid- 
ification by the technique of matched asymptotic expan- 
sions. This procedure has been detailed in several pub- 
lications for models that are similar to ours, and hence 
we will only outline the results. Each solid has to reject 
its minority component and the ternary impurity into 
the liquid in order to grow. Since the concentrations are 
locally conserved quantities, mass balance at the inter- 
face implies that they obey boundary conditions of Stefan 
type at the moving boundary, i.e. 



Dd n u = v n (m -u s ), 
-Dd n c = v n (l - K)ci. 



(27) 
(28) 



Here, v n and d n are the normal velocity of the interface 
and the derivative normal to the interface, and ui, u g , and 
ci are the values of the concentrations at the liquid and 
solid sides of the interface, respectively. These equations 
are valid for both solid-liquid interfaces; note that on 
the a-liquid interface, m — u s > 0, whereas on the (3- 
liquid interface, ui — u s < 0. The concentrations at the 
interface are related to temperature, shape, and speed of 
the interface by a generalized Gibbs-Thomson condition, 



T = T E T m u ui - rhci — DC — v n /fi k 



(29) 



where the upper (lower) sign is for the a [0) phase, the 
liquidus slopes m u and fh are given by Eqs. ( |l5| ) and 
(Jig), r is the Gibbs-Thomson constant, K. is the local 
curvature of the interface, /j,k is the linear kinetic coeffi- 
cient, and the concentrations on the solid side arc linked 
to those on the liquid side via Eqs. ([l0|) and (|l^). With- 
out the last (kinetic) term, Eq. (|29[) is a statement of 
local equilibrium at the interface, including capillary ef- 
fects. The Gibbs-Thomson constant T is given in physical 
units by r = j s iTe/L, where j s i is the surface tension 
of the solid-liquid interface, and L is the latent heat of 
melting. In the units of our model, this becomes 



j s iT E /a. 



(30) 
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The surface tension 7^ is obtained as in Ref. |30| by 
first solving numerically the one-dimensional stationary 
versions of the equations of motion to obtain the interface 
profiles, and then computing the excess free energy per 
unit surface by inserting the profiles into the free energy 
functional. Note that, since the free energy density is 
taken dimcnsionless here, surface tensions have units of 
length. The solid-solid surface tension can be calculated 
analytically since along the a/3-interface = 1, and we 
obtain 7ss = (2/3) W u . 

From the surface tensions, we can determine the con- 
tact angle 6, that is, the angle between the horizontal and 
the solid-liquid interfaces at a trijunction point where the 
solid-solid interface is vertical. Using Young's condition 
of mechanical equilibrium, we get 



sin 9 



Iss 



(31) 



We have not explicitly calculated the value of the 
kinetic coefficient /ik that appears in the last term of 
Eq. (p9|). This would require to compute several in- 
tegrals in the coupled variables u and <p through the 
solid-liquid interface (see Ref. Q for more details in a 
simpler case), which can only be done numerically. Fur- 
thermore, we have neglected in the above analysis other 
non-equilibrium effects, and in particular solute trapping 
|Q that is genera lly pre sent in phase-field models for 
alloy solidification 29 41J1 ■ It is known that solute trap- 
ping modifies the compositions on both the liquid and 
the solid sides of a moving interface. This generates cor- 
rection terms both in the Gibbs-Thomson condition and 
the mass balance relations, Eqs. (E8|) and (|2^). How- 
ever, these corrections are proportional to the interface 
velocity, and are expected to be small for the range of 
solidification speeds used in our present simulations. In- 
deed, as we will see below, the non-equilibrium effects are 
not entirely negligible; however, they are not important 
enough to justify a detailed analysis that would be quite 
involved 



III. LAMELLAR STEADY STATES 

We chose as a testing ground for our model the sim- 
ulation of lamellar steady-state solutions. This has the 
additional benefit of providing us with the initial config- 
urations needed for the simulations of large-scale arrays 
described below. In the laboratory frame, the sample is 
pulled with velocity v p in a constant temperature gradi- 
ent G along the z-axis. This means that in the sample 
frame, the isotherms move towards the positive z direc- 
tion with velocity v p . Consequently, the temperature at 
a given point (x, z) of the sample is 



T(z, t) = T E + Gz ~ v p t, 



(32) 



The equations of motion were simulated by an explicit 
Euler algorithm with timestep At on a simple square grid 
of spacing Ax using standard finite-difference formulae. 
For simplicity, we chose W u = = W. In the following, 
unless otherwise stated, lengths will be measured in units 
of W, time in units of r, and temperatures in units of Te- 
We chose D = D = l, a = l, a = l/8 (parallel eutectic 
solidus and liquidus lines), G — 0.001, v p between 0.005 
and 0.02, and various values of b and K, with a = and 
e s = —blnK. Since the equation for the composition u 
is of fourth order, the critical timestep for the occurrence 
of numerical instabilities scales as Ax 4 . The allowed grid 
spacing Ax, however, is limited by the requirement that 
the smooth interfaces be sufficiently well resolved to avoid 
strong numerical anisotropies and lattice pinning. We 
found that Ax = 1 and At — 0.025 provided a good 
compromise between efficiency and accuracy. 

The simulations were started with a single pair of flat 
lamellae in contact with the liquid in a box of lateral size 
A. The concentrations were set to the equilibrium values 
in each phase. For the subsequent evolution, periodic 
boundary conditions were used in the direction parallel 
to the temperature gradient, while the concentrations in 
the liquid were kept at fixed values Uoo and Coo at the up- 
per end of the simulation box. At the lower (solid) end, 
no boundary conditions are needed since the fields do 
not evolve. During the runs, the simulation box was pe- 
riodically shifted to follow the interface. Convergence to 
the steady-state solution was checked by computing the 
average change of the phase field in the moving frame 
during the advance of the isotherms by one lattice spac- 
ing. Furthermore, the interface shapes (given by the level 
set <fi — for the solid-liquid interface and by u = in 
the solid, that is for <f> > 0, for the solid-solid interfaces) 
are extracted by interpolation of the fields between the 
lattice points. This procedure yields a resolution far su- 
perior to the grid spacing. The average undercooling of 
the interface is then 



AT{t) = T E - T int {t) = -G 



1 



C(x, t)dx — v p t 



where we have chosen the origin of the z-axis at the eu- 
tectic isotherm for t = 0. 



(33) 

where £(x, t) is the z-position of the extracted solid-liquid 
interface as a function of x at time t. The simulations 
were stopped when the undercooling was to within 10~ 4 
of its extrapolated final value. 

We first discuss the special case of a pure (binary) eu- 
tectic at the eutectic composition, — Cqq — (note 
that we omit the impurity terms in the free energy and 
the equation of motion for the impurities when Coo = 0). 
For our symmetric phase diagram, there is no global dif- 
fusion boundary layer in this case, and the diffusion field 
in the liquid decays exponentially on a scale of A. Hence, 
a box length parallel to the temperature gradient of about 
five times A was sufficient to obtain results that are inde- 
pendent of the box size. The interface relaxes exponen- 
tially to its steady state, with relaxation times of order 



G 



A 2 / D; on a typical modern workstation, the convergence 
takes a few hours. 

In contrast, for Uoo ^ and/or Coo 7^ 0, solute redistri- 
bution leads to a boundary layer of thickness Id — D/v p , 
much larger than A. Hence, box sizes along the growth 
direction of several times Id have to be used; in addi- 
tion, the interface position now follows a damped os- 
cillation with exponential envelope and decay times of 
the order D/v 2 , much larger than X 2 /D. As a result, 
when a boundary layer is present the convergence of a 
run takes several days of CPU time. The oscillatory re- 
laxation of the interface position is compatible with the 
Mullins-Sekerka dispersion relation at zero wave vector 
that predicts a complex decay rate. 

Let us compare our results to the well-known Jackson- 
Hunt (JH) relation between lamellar spacing and inter- 
face undercooling j| , generalized to include the effect of 
the ternary impurities, 



AT(A) = ATjh{\) + rhcoo/K, 



A-Tjh — -AT m i n 



A A. 



A,. 



A 



(34) 



(35) 



The curve AT versus A exhibits a minimum at a spacing 
Amin, where 



AT • — 



2m u Au Tsm9P(r])i 



77(1 — rj) y 2Dm u Au ' 



Amin — 



/ 2rsin6»T> 

m u Auv p P(rj) 



(36) 



(37) 



Here, Aw = up — u a = 2 is the concentration differ- 
ence between the two solids at the eutectic temperature, 
?/ = (iioo — u a ) / Au is the volume fraction of /3-phase in 
the solid, P(r]) = Y^=i sin 2 (nnn) / (mi) 3 , and 9 is the 
contact angle defined by Eq. (|3l|). For W u = W$ — 1, 
we obtain numerically y s l = 1.04, which together with 
7 SS =2/3 gives » 19°. Note that Eqs. @ and (§7|) 
are valid only for our choice of a symmetric phase dia- 
gram; see Ref. |23| for a discussion of the general case. 
It should be kept in mind that the JH theory is approx- 
imate since it uses a flat interface to calculate the diffu- 
sion field. Nevertheless, it has been shown by boundary 
integral simulations that the error is small for small 
contact angles 6 and close to the spacing A m i n , such that 
it can be used as a semi-quantitative test for our phase- 
field model. 

We computed the interface undercooling in our model 
for various pulling speeds and two different values of the 
mobility exponent n in Eqs. (|2^) of the mobility func- 
tions. Let us first discuss the results for n = 1, which 
corresponds to the simplest form of the mobility that 
has been widely used before. The simulated undercool- 
ings are slightly higher than the JH prediction, but the 
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FIG. 2. Average interfacial undercooling versus lamellar 
spacing for several values of the pulling speed v p and the 
mobility exponent n. Lines: prediction of the Jackson-Hunt 
theory, Eq. (^); symbols: simulation results. Filled circles 
correspond to steady states that are unstable with respect to 
1-A oscillations. 



overall shape of the curve is perfectly reproduced. The 
difference can be attributed to the non-equilibrium ef- 
fects (interface kinetics, solute trapping) present in the 
phase-field model, but neglected in the JH theory. In- 
deed, the differences between our simulations and the 
JH prediction are larger for higher v p . Furthermore, we 
obtained A m in and AT m i n by fitting our simulation re- 
sults to Eq. (|35|) and found that the scaling relation 
A min u p = const, that can be derived from Eq. (|37| ) is 
well satisfied. Regarding the impurity contribution to 
Eq. (|34]), we conducted simulations for various impurity 
concentrations, impurity liquidus slopes and partition co- 
efficients and found good agreement with the predicted 
behavior. In particular, we verified that the spacing A m i n 
was not appreciably modified by the addition of impuri- 
ties. 

The range of lamellar spacings that can be simulated 
is limited by two effects that are intrinsic to our model. 
For spacings smaller than ~ 16W (~ 8W for each in- 
dividual lamella), the diffuse interfaces at the two sides 
of a lamella start to overlap, which leads to strong cor- 
rections to the sharp-interface limit and ultimately to 
lamella elimination. 

For too large spacings, in turn, new lamellae of the op- 
posite phase nucleate in the centers of the initial lamellae, 
leading to a lamellar array with one third of the initial 
spacing. This is the result of a "spinodal decomposition" 
that takes place in the interface. Indeed, the equation 
for the composition in the solid far from the interface is 
exactly the classical Cahn-Hilliard equation |42[] , which 
is known to exhibit phase separation without nucleation 
in a composition range where the free energy density has 
a negative curvature (d 2 f so \/dc 2 < 0). Far inside the 
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solid, this has no importance here because the mobility 
is zero and hence no dynamics takes place. Well within 
the liquid, there is no unstable concentration range since 
the liquid free energy has a single-well structure. But in 
the diffuse interface, new domains may form when the 
concentration falls within the unstable range. Accord- 
ing to the JH theory, the deviations of the concentration 
from the equilibrium value at the interface scale as Pe A, 
where the Peclet number Pe = X/Id = Xv p /D; hence, 
the maximum spacing A max that can be simulated before 
nucleation sets in increases as v p decreases. Indeed, we 
find A max /W ~ 28 for v p = 0.02 and X max /W ~ 58 for 
v p =0.01. 

It seems useful at this point to comment on the impli- 
cations of these limitations for the choice of the computa- 
tional parameters for large-scale simulations. The range 
of initial lamellar spacings of interest for the present 
study ranges from A m in to about 1.5 A m i n . Since we 
want to simulate the linear instability of a lamellar front, 
which may involve considerable variations of the local 
lamellar spacing, the model should work for a sizeable 
range of spacings, say at least for spacings that are ±50% 
of the initial value. This means that we must have 
X m in/W > 32 because of the low-spacing limitation. 

Next, A m in should not be much larger than this value, 
since the computer time necessary to simulate the evolu- 
tion of an array of initial spacing A m ;„ can be estimated to 
scale as A min (number of grid points: A m i n x A m i n ; time 
for the interface to advance by one spacing: A m j n /u p ; 
using A min u„ =const., we get t C pu ~ A min ~ v p ~ 5/2 ). 
From Fig. @, we can see that n — 1 and v p = 0.01 
give A m i n of the right order of magnitude; however, since 
Amax/W^ ~ 58, the available range of lamellar spacings 
is somewhat small (for an initial spacing of 1.5 A m j n , nu- 
cleation would set in for an increase of the local spacing 
by only 30%). Since the range of available spacings in- 
creases with decreasing pulling speed, one possible solu- 
tion would be to further reduce v p . However, as discussed 
above the necessary computer time rapidly becomes pro- 
hibitive. 

Another way out is to change the exponent in the equa- 
tions for the mobilities, Eqs. (p2|). If we choose n > 1, 
the diffusivity is increased in the whole interfacial re- 
gion, whereas it remains zero in the solid. This leads to 
higher diffusion currents along the surface than for n = 1 . 
Hence, the pileup of the rejected atoms at the interface 
is lower, and consequently A max is higher. The price to 
pay is that this model, in its sharp-interface limit, is not 
equivalent to the classical JH model, but contains addi- 
tional surface diffusion terms ^3|. However, as shown in 
Fig. ||, the qualitative behavior of the undercooling ver- 
sus spacing curve does not change. For v p = 0.01, A m ; n 
is larger than the theoretical JH-value by about 10%, 
whereas AT m j n is about 15% too low. On the other hand, 
Amax/W^ = 74, such that we now have at our disposal a 
sufficient range of spacings. 

For these parameters, we observe for large spacings the 



well-known period-preserving oscillatory instability that 
sets in at about 2A m j n [Q. Even beyond the threshold of 
this instability, steady states can be reached to within an 
excellent precision, because we start from an exactly sym- 
metric initial condition and because the numerical noise 
of the phase-field approach is extremely low. To trigger 
the instability within a reasonable simulation time, an 
explicit perturbation that breaks the symmetry between 
the two phases had to be added. Such unstable steady 
states are shown as filled symbols in Fig. |[ 

The mechanism for lamella creation by nucleation is 
in fact very useful for the simulations of well-developed 
colonies where lamellae are frequently created at the so- 
lidification front. We want to confront our simulations to 
the experimental findings of Akamatsu and Faivre jl8) , 
who work with thin samples of a transparent eutectic 
alloy enclosed between parallel glass plates. In their ex- 
periments, creation of new lamellae takes indeed place 
predominantly in the center of already existing lamel- 
lae. However, the detailed mechanism is still unknown. 
New lamella do not form by nucleation, since the inter- 
facial undercoolings are not high enough. Most likely, 
the "pockets" in the center of large lamellae are "in- 
vaded" from pre-existing neigboring lamellae of the op- 
posite solid by fingers that grow in the meniscus between 
the glass plates and the growing solid. The point here 
is that the modeling of such a process is out of reach for 
our present computational resources, since it is neces- 
sarily three-dimensional. Within the framework of two- 
dimensional simulations, we simply need a criterion to de- 
cide when new lamella should form, and the "automatic" 
implementation of a maximal lamellar width A max in our 
model is an adequate solution that avoids the implemen- 
tation of an explicit nucleation rule, as done for example 
in Rcf. @. 

IV. LINEAR STABILITY OF LAMELLAR 
ARRAYS 

A. Theory 

We have recently performed a detailed linear stability 
analysis of a lamellar eutectic interface in the presence 
of ternary impurities. Rather than to repeat the calcu- 
lations here, we will give a brief summary of the main 
assumptions and results before discussing the phase-field 
simulations. Our analysis is an extension of the method 
used by Datye and Langer (DL) to analyze the stability 
of lamellar arrays without impurities [^J . It is based on a 
perturbation scheme for the Jackson-Hunt solution and 
proceeds as follows. 

1. The coordinates of the trijunction points are chosen 
as fundamental variables to describe the state of the 
perturbed system. This amounts to a "discretiza- 
tion" of the original continuous system. Each tri- 
junction point has two degrees of freedom, namely 



8 



its x and z positions (motion parallel and normal 
to the isotherms, respectively). 

2. For a lamellar interface that is gently curved on a 
scale much larger than A, the lamellae are assumed 
to grow perpendicular to the envelope of the com- 
posite front (Cahn's hypothesis). This connects the 
time derivative of the local lamellar spacing to the 
shape of the front. For example, in a protrusion 
where the front curves outward, the local spacing 
increases during further growth. 

3. Given the positions of the trijunction points, the ac- 
tual interface shape is replaced by a piecewise pla- 
nar interface, and a perturbed diffusion field is cal- 
culated. The Gibbs-Thomson equation is then used 
to obtain an eigenvalue problem for normal modes, 
i.e. perturbations proportional to exp(ifca; + cut), 
where k is the wave vector of the periodic per- 
turbation, and w its growth rate. The solutions 
of the eigenvalue equation give the dispersion rela- 
tions u>(k). Since there are four degrees of freedom 
per lamella pair (two for each trijunction), u>{k) 
has four branches. Of those, there are two that are 
relevant for the long-wavelength instability we are 
interested in. 

It turns out that the final result of this rather compli- 
cated analysis can be understood in terms of an effective 
front approach. Namely, one can separate two scales: the 
local lamellar spacing, and the large-scale smooth enve- 
lope of the lamellar front. The evolution of the local spac- 
ing is slaved to the shape of the front by the assumption 
of normal motion (Cahn's hypothesis). On the scale of 
the smooth front, the lamellar structure introduces an in- 
terfacial undercooling that is approximately given by the 
JH law, taken with the local spacing and interface veloc- 
ity. Using these ingredients, it is possible to include the 
lamellar geometry in the usual Mullins-Sekerka stability 
analysis and to obtain the dispersion relation. This re- 
sult can be recovered from the more complicated discrete 
analysis, with one additional ingredient. The eutectic dif- 
fusion field that governs the exchange of atoms between 
neighboring lamellae gives, when perturbed on a large 
scale, a stabilizing contribution to the total interfacial 
undercooling. The functional form of this stabilization 
is the same as for the surface tension terms, and this ef- 
fect can therefore be included in the simple effective front 
approach by simply "renormalizing" the capillary length. 

The two main results of this analysis are that (i) the 
instability threshold is close to the well-known constitu- 
tional supercooling criterion, with a small capillary cor- 
rection, and (ii) in contrast to the Mullins-Sekerka in- 
stability, where unstable modes always have real growth 
rates, the lamellar structure may lead to complex growth 
rates, and hence to oscillatory modes. The origin of these 
oscillations can be understood as follows: in a protrusion 
of the front, the lamellar spacing increases. This leads to 
a local change in the JH undercooling that, for a small 



distortion of an array of spacing Ao, is proportional to 
the slope of the JH plot. For A > A r „in, this provides 
a "restoring force" for the large-scale front. Since only 
the change with time of the lamellar spacing (but not the 
spacing itself) is related to the shape of the front, the 
dispersion relation becomes quadratic in lu, instead of 
the linear Mullins-Sekerka equation. There are two solu- 
tions to this equation for each wave vector k. In physical 
terms, this is the consequence of the additional "internal 
degree of freedom" A of the front. As discussed in detail 
in Ref. p3| , real and complex growth rates may occur, 
depending on the ratio G/v p , the lamellar spacing, and 
the impurity content. For large enough spacings, when 
the "restoring force" mentioned above is strong enough, 
we expect that the complete dispersion relation is com- 
plex. One of the goals of the present paper is to test this 
prediction by direct simulation of the basic equations of 
motion. 



B. Single mode simulations 

Let us first study the behavior of a single unstable 
mode of a lamellar array with initial spacing Ao- The 
parameters besides Ao that control stability are the im- 
purity content and the ratio G/v p . We define the dimen- 
sionless parameters 



w = 



A = A /A n 



mAc mCooll/K — 1) 



m u Au 



2DG 
v p m u Au 



mAc 



2DG 
v v mAc 



(38) 



(39) 



(40) 



Here, A m ; n is obtained from an interpolation of our sim- 
ulation data shown in Fig. [| For n = 4 and v p = 0.01, 
Amin ~ 34. The freezing ranges of the eutectic and 
the impurities are, expressed in the parameters of our 
model, m u Au = 2Te/u, and mAc = m(l/K — l)coo = 
6(l-X) 2 T B c 0O /(ifa). 

A lamellar array is prepared by replicating the steady- 
state solution for one lamella pair N times. We apply a 
cosine perturbation to the steady state and impose the 
initial condition 



(x, z, 0) = 4>o(x, z + A cos(2nx/N)), 



(41) 



where 4>q{x,z) is the steady-state solution. The other 
fields (u and c) are perturbed in the same manner. The 
perturbation amplitudes Aq are usually much smaller 
than the interface width (typically, Aq/W < 0.1), and 
the values on the grid points are obtained by linear in- 
terpolation of the numerical steady-state solution. 

To analyze the evolution of the system, we store peri- 
odically the positions of all the interfaces (solid-solid and 
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FIG. 3. Evolution of oscillatory modes. Growth direction 
is from bottom to top, and three successive frames are shown 
from left to right. Shown are the solid-solid-interfaces, as 
well as successive snapshot pictures of the solid-liquid inter- 
faces. The system is perturbed with a single cosine mode of 
dimensionless wave vector k = 0.2. Simulation parameters: 
v„ = 0.01, G = 0.0005, A = 40, K = 0.5, c x = 0.08, b = 10, 
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solid- liquid) . In addition, we determine the positions of 
all the trijunction points by searching the intersections of 
the level curves = and u = 0. The coordinates of the 
trijunction point to the left of the ^-lamella (y = a, (3) in 
the lamella pair number i are labeled [x\,z"). We define 
the deviations of the trijunction point coordinates from 
their steady-state values, 



y, 



z 

•JC? 



as well as their discrete Fourier transforms, 
1 N 

X v {n, t) = — ^2 exp(27riKn) 

n=l 
1 N 

Y v (n,t) = — ^y^exp(27riKn), 



(42a) 
(42b) 



(43) 
(44) 



n=l 



where k — fcAo/(27r) is a dimensionless wave vector. 

In Fig. [|, we show the evolution of an array of five 
lamellae, started in a single mode with n — 0.2. An 
oscillatory 5-A mode develops. Its amplitude grows un- 
til a lamella termination occurs at z/\$ — 50. Subse- 
quently, the system shows a decaying 4-A-oscillation and 
approaches a steady-state solution with 4 lamella pairs. 
To analyze this evolution, we use the Fourier components 
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FIG. 4. Real part of the Fourier amplitude X a (K,t) as a 
function of time, for the run of Fig. |^. Also shown are the 
fits to growing (decaying) oscillating exponentials. 

X v {n,t). Since the initial perturbation is not propor- 
tional to the (unknown) eigenvector corresponding to a 
single mode of the complete (continuous) system, the 
Fourier components for n ^ 0.2 will not remain zero. 
However, they remain sufficciently small to be neglected 
in the data analysis. In Fig. we show the evolution of 
Ke[X a (K,t)] versus time, for k — 0.2 before the lamella 
elimination, and for k = 0.25 afterwards (note that the 
elimination of one lamella pair corresponds to a change in 
the unperturbed lamellar spacing Ao). Oscillating modes 
correspond to complex growth rates. We define the di- 
mensionless growth rate by 



(45) 



with Q r and fl; real. The growth rate is determined by 
a fit of the data to the function 



Re[X a (K,t)} = Aexp(ft r i) sin[Oi(f - t )], 



(46) 



where t is measured in units of \q/v p . In practice, we ob- 
tain a value of t Q , i.e. the time of one of the zero crossings, 
by numerical interpolation, and then use a least-squares 
fitting procedure with A, Q r , and £li as free parameters. 
As can be seen from Fig. |[ the fit is excellent. Sur- 
prisingly, the fit remains accurate up to the immediate 
vicinity of the lamella termination event. This indicates 
that the system is well described by a single, exponen- 
tially growing mode even for large deformations of the 
initial array. In particular, the linearization that is the 
basis for the theoretical analysis remains valid even if the 
lateral displacements are large, i.e., y\ /Aq ~ 1. 



C. Dispersion relations 

The simulation and fitting procedures outlined above 
were carried out for various values of the control param- 
eters and arrays of different sizes to construct the dis- 
persion relations J7(k). In Fig. we show a comparison 



10 



a 



0.3 



0.2 



0.1 



0.0 - 




-0.1 



0.0 



0.1 0.2 

K 



0.3 



a 



0.2 



0.1 



0.0 



(b) 



O 



G 0.2 - ^j" 

0.0 



O 



0.1 0.2 
K 



-0.1 



0.0 



0.1 0.2 

K 



0.3 



FIG. 5. Plots of the dimensionless growth rate f2 versus 
dimensionless wave vector k. The main graphs show the real 
part Q r , the insets the imaginary part Qi for (a) A = 1.175, 
w = 0.2, g = 0.05 and (b) A = 1.175, w = 0.2, g = 0.1. Filled 
symbols and lines: theoretical predictions from Ref. pj| ; cir- 
cles: real modes (O, = 0), squares: complex modes (Qi 7^ 0). 
Open symbols: simulation data; triangles: real modes, dia- 
monds: complex modes. 



of the obtained data with the theoretical predictions of 
Rcf. Q for two different values of the temperature gra- 
dient. For both dispersion relations, there are stationary 
(fi real) and oscillatory (f2 complex) modes. According 
to theory, for g = 0.05 the fastest growing mode is sta- 
tionary, whereas for g = 0.1 is is oscillatory. 

In all cases, the nature of the mode (stationary or os- 
cillatory) agrees with the theoretical predictions. Fur- 
thermore, the oscillation frequency of the complex modes 
(f2i) is always in good quantitative agreement with the- 
ory. In contrast, the growth rates (f2 r ) are in good agree- 
ment only for small wave numbers; for large wave num- 
bers, the simulated growth rates are systematically much 
smaller than predicted by theory, and the difference in- 
creases with the dimensionless wave vector. Therefore, 
in the simulations at g = 0.1, the fastest growing mode 
is stationary, and not complex as predicted by theory. 
For A = 1.47 and g — 0.1, we obtain a stability spectrum 
that is entirely complex (data not shown), both in theory 
and simulations. 

Just as the JH theory, our stability analysis of a lamel- 
lar array contains several simplifying assumptions. It 
is therefore necessary to check whether the differences 
between theory and simulations are due to the approx- 
imations made in the stability analysis, or due to the 
phase-field approach, which is a genuine representation 
of the original free-boundary problem only in the limit 
W/\ — » 0. Therefore, we focused on a single complex 
mode at g = 0.1 and k — 0.2 (5-A-oscillation) and con- 
ducted a series of runs with decreasing pulling speed v p . 



Since A n 



-1/2 



, we increased the spacing Aq to keep 



the reduced spacing A constant. The temperature gradi- 
ent G was also decreased to keep g constant. 

The results for the growth rate fi r versus v p are shown 
in Fig. The data fall on a straight line, and by extrap- 
olation to v p = we find fl r (v p = 0) = 0.085. In con- 
trast, the variation of the oscillation frequency is very 
small (from = 0.291 at v p = 0.01 to Q; = 0.302 at 
v p = 0.005). This linear variation of Q r with v p indi- 
cates that the dominant corrections to the sharp-interface 
limit of the phase-field model scale as W/Id — Wv p /D. 
Corrections in the other involved small ratios, W/\q 
and Xq/Id, seem to be subdominant, since both scale as 
1/tJVp~ at constant A. An example for a correction that 
scales as W/Id is the interface kinetics; however, insert- 
ing a kinetic term in the Mullins-Sekerka analysis does 
not lead to a linear variation of the growth rate with the 
kinetic coefficient. The solute trapping effect also scales 
as W/Id, but since it is quite involved to evaluate its 
influence on the growth rates, we have not investigated 
this issue in more detail. We checked, however, that the 
variation of f2 r with v p is not a consequence of the surface 
diffusion term introduced by our choice n = 4 in the mo- 
bility function: a simulation with n = 1 and comparable 
A yielded similar results. 

The simulated growth rate, extrapolated to v p = 0, 
is still markedly different from the theoretical prediction 
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FIG. 6. Real part of the dimensionless growth rate Q versus 
pulling speed for A = 1.175, g = 0.1, w — 0.2, and k — 0.2. 
Symbols: simulation result. Dashed line: linear extrapolation 
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FIG. 7. Sideways velocity, dy/dt, of one trijunction point 
for A = 1.175, g = 0.1, w = 0.2, « = 0.2, and v p = 0.00751. 
Solid line: data extracted from the simulated curve y(t). Dot- 
ted line: prediction of Cahn's hypothesis. Dashed line: best 
fit to Eq. (0). 



Q r = 0.1365. We therefore checked several assumptions 
that are used in the linear stability analysis, in particular 
Cahn's hypothesis that the lamellae always grow perpen- 
dicular to the large-scale front. Expressed in terms of the 
trijuction point coordinates defined in Eq. ([42]), this as- 
sumption reads 



dtyt 



.^p. e<* 



in) 



(47) 



for the trijunction point to the left of the nth (3 lamella. 
From the simulation data, explicit values oidty^ and the 
vertical displacements are available, and Eq. j47| ) 

can be directly checked. As shown in Fig. [7, Cahn's hy- 
pothesis is clearly violated. We tried to fit the difference 
of the right-hand-side and left-hand-side of Eq. (E^) to 
various functions of the trijunction coordinates and found 
that the modified equation 



(48) 



yields a good fit with a single adjustable parameter J5, 
as is shown in Fig. |7[ 

We repeated the above fit for all our simulation data. 
Rather remarkably, the correction given by Eq. ( f4g| ) 
works both for oscillatory and stationary modes, and the 
fit parameter B behaves smoothly in the crossover re- 
gions. This shows that the violation of Cahn's hypothesis 
is a consequence of the local front geometry, and not a 
cooperative effect depending on the nature of the mode. 
We found that the fit parameter B mainly depends on the 
wave vector k and the reduced spacing A. We also found 
weak dependencies on the temperature gradient g and 
the impurity partition coefficient K. More importantly, 
B is almost independent of the ratios Ao/W and W/Id- 
for the series of runs of Fig. ||, B varied between 0.240 
for v p = 0.01 and 0.215 for v p = 0.005. Even though B 



decreases with v p , it does not extrapolate to 0, such that 
this effect is not an artefact of the phase-field model. In 
Fig. H we show the fitted values of £?, rescaled by the 
reduced spacing A, versus the wave number k for various 
external parameters. Neglecting the weak dependencies 
on K and g, a fit to all the data points yields 



B = B Ak 2 



(49) 



with B = 4.96. 

Although Cahn's rule is violated, the resulting devi- 
ations of the growth angles from 90° are very small. 
To see this, let us use the geometrical relation d t y^ = 
— tan <5^, where 5@ is the angle between the solid-solid 
interface at the trijunction and the z direction. From Eq. 
(|48| ) we can calculate the deviation of 6@ from the value 
predicted by Cahn's rule, which is (£"+1 — £")/^o- I n 
our simulations, this deviation never exceeded 1°. Due 
to the finite interface width of the phase-field model, it 
is very difficult to measure angles directly at the trijunc- 
tion points, and such small deviations cannot be resolved. 
Therefore, the procedure outlined above that uses the 
whole trajectory of a trijunction point is the only way 
to obtain quantitative information about the violation of 
Cahn's rule directly from the simulations. It should be 
emphasized, however, that while the deviation itself is 
small, since the growth angle 5@ is itself small, the ra- 
tio of the two is not necessarily small. Indeed, it can 
be seen from Fig. (?] that the correction constitutes a 
sizeable fraction of the growth angle. This explains why 
such a small deviation can induce quite large shifts in the 
stability spectrum. 

The functional form of Eq. ( f49| ) allows us to draw sev- 
eral interesting conclusions. First, since the coefficient 
B is proportional to the reduced eutectic spacing A, but 
almost independent of the impurity content, this effect 
is not specific to ternary alloys, but should also occur in 
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FIG. 8. Fit parameter B divided by the reduced spacing A 
versus wave vector k for various sets of control parameters. 
The line is a best fit of all data with Eq. 



binary eutectics. Secondly, a deviation from Cahn's rule 
has been previously reported in binary eutectics J44] , [45|] . 
However, this deviation becomes important only for very 
strong temperature gradients, i.e. for j » 1, when the 
front is almost flat, except in the immediate vicinity 
of the trijunction points. Here, we are in the opposite 
regime, with g <C 1, and a front that consists of round 
arcs between the trijunction points; since, in addition, B 
is almost independent of G, we conclude that the effect 
described by Eqs. ( fig ) and fl49| ) seems not to be captured 
by these calculations. 

Additional physical insight can be gained by inserting 
back Eq. ( [49| ) into Eq. (ff8|). The correction to Cahn's 
rule is proportional to n 2 y. In a continuum limit, where 
the function y(x) is a smooth interpolation of the lateral 
trijunction displacements, this corresponds to a second 
derivative, d xx y(x). For y varying slowly on the scale of 
Ao, we have A(x) ~ Ao(l + d x y), such that the correc- 
tion to Cahn's rule is proportional to the gradient of the 
local spacing. The motion of the trijunctions is there- 
fore a combination of the perpendicular lamellar growth 
considered before and a small lateral drift proportional 
to the gradient of the local spacing. While such a term 
certainly appears to be reasonable, a more detailed the- 
oretical analysis is clearly warranted. 

This violation of Cahn's hypothesis explains the re- 
maining discrepancies between our simulation results and 
the theory. To modify the theory by the inclusion of the 
corrective term in Eq. (H) seems possible, but is out of 
the scope of the present paper. 



V. DYNAMICS OF COLONY FORMATION 

To study the instabilities that lead to the formation 
of colonies, we constructed large arrays as described be- 
fore, and perturbed the steady-state solution by a spa- 



tial displacement of the fields along the z direction. The 
amplitude of the displacement was a random variable of 
x with a white noise spectrum and an amplitude com- 
parable to one lattice spacing. The goal was to study 
the initial instability of such random arrays as well as 
the nonlinear dynamics of well-developed colonies. The 
latter required long runs in big systems. The necessary 
computational power was attained by porting our sim- 
ulation code on a parallel CRAY T3E computer. We 
used a simple domain-decomposition scheme for paral- 
lelization, i.e. every processor calculated a part of the 
system. A load-balancing algorithm that adjusted the 
domain boundaries as a function of the computational 
load for each processor was used to optimize the yield. 

For = (eutectic composition), the initial evo- 
lution of the lamellar array is a linear superposition of 
the long-wavelength modes desribed in the previous sec- 
tion. That is, if we decompose the set of trijunction dis- 
placements into Fourier modes, each mode grows with 
the (real or complex) growth rate that was determined 
in the single-mode simulations of the preceding section. 
In Fig. ^ we show the resulting evolution for the same 
control parameters as in Fig. |(a). The fastest growing 
mode is real with a wavelength of about 12Ao- Indeed, 
this mode dominates the interface shape in the second 
snapshot, where the first lamella termination events have 
occured. At later times, the linear description becomes 
evidently invalid. The further evolution is characterized 
by the growth of long protruding fingers, as can be seen 
in the last snapshot picture. These fingers, however, do 
not reach a steady-state configuration up to the end of 
our simulation: their shape continuously changes, and 
there are some tip-splitting and overgrowth events. To 
highlight this feature, we show in Fig. |l^ a complete 
plot of the whole solidified sample, where we have omit- 
ted the greyscale for clarity, and where we have marked 
the trajectories of the "deep grooves" between neighbor- 
ing fingers. This run was performed on a lattice of size 
1600 x 1200 and totals 15 x 10 6 iterations. On the CRAY 
T3E, this run required about 3000 hours single processor 
CPU time. 

In Fig. [0], we show a run with again 40 lamellae, but 
now with both a larger temperature gradient and a larger 
initial spacing. Under these conditions, the instability de- 
velops more slowly, and the dispersion relation is entirely 
complex, such that we expect propagating or oscillatory 
modes. Indeed, on the right side of Fig. |ll|, there is 
an oscillatory "breathing mode" with wavelength about 
10Ao, whereas on the left side, a travelling perturbation 
of the lamellar pattern can be seen. The run was not con- 
tinued after the first lamella termination events, since the 
nonlinear regime is expected to lead to similar fingered 
patterns as in Fig. |[ 

A quite different scenario occurs for off-eutectic com- 
positions. An example is shown in Fig. |l2|. The linear 
regime is still in good agreement with the predictions 
of Ref. 1 23 1 . In particular, for sufficiently off-eutectic 
compositions, the impurity-induced long-wavelength in- 
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FIG. 9. Snapshot pictures of a run with 40 lamella pairs at 
eutectic composition and g — 0.1, A = 1.175, and w = 0.2. 
From top to bottom: t/r = 0, 100000, 375000. In the solid, 
red and blue represent the two solid phases. In the liquid, the 
green intensity is proportional to the impurity concentration; 
the small blue and yellow "halos" in advance of the grow- 
ing lamellae are a visualization of the interlamellar (eutectic) 
diffusion field. 




FIG. 10. Global view of the same run as in Fig. y, without 
greyscale. Thin lines: solid-solid interfaces. Thick solid line: 
final solid-liquid interface. Thick dashed lines: trajectories of 
the grooves between fingers. There are two tip-splitting and 
one finger overgrowth event. Note the concave part of the 
final front in the center of the leftmost finger: a tip-splitting 
event will soon take place. 
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FIG. 11. Run with 40 lamella pairs at eutectic composition 
and g = 0.2, A = 1.5, and w — 0.2. The dispersion relation 
is entirely complex, and oscillatory patterns appear. Thin 
lines: solid-solid interfaces. Thick solid line: final solid-liquid 
interface. 



stability competes with the 2A-oscillatory instability that 
is already present in binary eutectics. For the temper- 
ature gradient and impurity content chosen in our ex- 
ample, the long-wavelength instability is stationary and 
faster than the 2A-0 instability. Indeed, we find that 
the Fourier spectrum of the trij unction displacements 
is initially dominated by the smooth long-wavelength 
modes, while the 2A-0 instability develops much more 
slowly. However, as soon as the instability becomes "vis- 
ible", that is, the amplitude of the perturbation exceeds 
~ 0.1A , localized finger-like structures develop around a 
lamella of the minority phase and rapidly grow ahead of 
the front. The fine lamellae act almost as "guides" for the 
well-developed fingers during the subsequent evolution. 
In particular, note the long minority lamella that is like 
a "spine" for the rightmost finger in the third snapshot 
(we remind the reader that we use periodic boundary con- 
ditions in the lateral directions; hence, this is not a "wall 
effect"). These structures, however, are only transient. 
In the final stage, when the colonies are well-developed, 
they have rather flat tops and sharper "corners" than 
the fingers at eutectic composition. In the flat parts at 
the center of the colonies, sometimes a period-doubling 
oscillatory mode develops until it generates some new 
lamellae and dies out. 

Structures such as the initial localized fingers are evi- 
dently non-linear. It thus appears that the linear regime 
of the instability is much shorter for off-eutectic than for 
eutectic compositions. It is presently unclear what pre- 
cisely triggers the formation of such fingers, and under 
which conditions they can form. In view of the necessary 
computer time, we did not carry out a detailed study to 
clarify these issues. 

VI. CONCLUSIONS 

We have presented a phase-field model for eutectic so- 
lidification in the presence of ternary impurities. This 
model has enabled us to carry out large-scale simula- 
tions of colony formation starting from arrays of up to 
40 lamellae pairs. 

In the linear regime, i.e. for small perturbations of 
the unstable steady-state growth front, these simulations 
have allowed us to critically test our previous linear sta- 
bility analysis p3fl . We find a good overall agreement 
with our theoretical predictions. Furthermore, a detailed 
treatment of the simulation data has allowed us to check 
the assumptions made in the linear stability analysis, and 
to precisely pinpoint the reasons for the differences be- 
tween the theory and simulation results. 

The most interesting conclusion is that the growth of 
the lamellae is not exactly normal to the large-scale enve- 
lope of the composite interface, a rule originally proposed 
by Cahn and used in the subsequent stability studies by 
Datye-Langer || and ourselves jg3) . This effect seems to 
be qualitatively different from the corrections to Calm's 
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FIG. 12. Snapshot pictures of a run with 20 lamella 
pairs at off-eutectic composition = 0.3), g = 0.1, 

w — 0.2, and Xq/W = 56. From top to bottom: 
t = 85000, 105000, 125000, 225000 (in units of r). 



rule reported previously in other theoretical studies of 
binary eutectics submitted to a strong temperature gra- 
dient Q,^]. The motion of the trijunction points can be 
roughly understood as a superposition of normal motion 
as stipulated by Cahn's rule and a slow "sliding" of the 
trijunctions along the front with a sideways velocity that 
is proportional to the gradient of the local lamellar spac- 
ing. The resulting deviations of the growth angles from 
90° are very small (below 1°); hence, a direct measure- 
ment of this effect in experiments is impossible, since a 
precise measurement of the growth angles is complicated 
by crystallographic effects, in particular the anisotropics 
of the surface tensions |lEj ]. However, the growth rates of 
the long-wavelength modes are very sensitive to a small 
change in this angle. Studying such modes can there- 
fore offer the possibility to experimentally test our re- 
sults. In particular, the consequences of this effect for 
the long-wavelength instability of binary eutectics will 
be discussed in a forthcoming study. 

Regarding the dynamics of fully developed colonies, 
we find that after the destabilization of the planar front, 
the array of two-phase cells undergoes a complicated and 
seemingly chaotic sequence of tip-splitting and cell elim- 
ination events. We were unable in our simulations to 
attain a steady-state configuration of the large-scale pat- 
tern, that is, the envelope of the front. This result is 
consistent with the fact that monophase cellular arrays 
in directional solidification of dilute alloys are unstable in 
the absence of crystalline anisotropy J37]|36|]. In fact, the 
lack of stability of the eutectic colonics in the absence of 
anisotropy suggests that the large-scale composite eutec- 
tic interface behaves qualitatively as a monophase front 
even beyond the linear regime. In this analogy, the ad- 
dition of solid-liquid or solid-solid anisotropy could po- 
tentially produce an effective anisotropy of the compos- 
ite interface that stabilizes its large-scale envelope. The 
quantitative exploration of this analogy, however, is far 
beyond the scope of the present work. 

Regarding the comparison between our simulations 
and the experimental observations of Ref. |l8[| , we find 
many similarities. In particular, we find in the simula- 
tions the oscillatory unstable modes predicted by our sta- 
bility analysis. Such wavy structures are also observed in 
the experiments. We also find that well-developed two- 
phase cells do not seem to reach a steady-state up to the 
largest times simulated. This is in agreement with the ex- 
periments, where no steady state has been reached even 
on length and time scales far superior to the range of our 
simulations (compare our Fig. [l(] to Fig. 14 of Ref. jlq]). 

A number of experimental observations, however, re- 
main to be understood. Firstly, unstable modes in the 
experiments are sometimes manifested as waves that are 
emitted by localized perturbations, such as grain bound- 
aries. These waves can propagate along the front, which 
remains planar, rather than be a transient that pre- 
cedes colony formation. Some of these propagating waves 
seem to have characteristics of solitary waves. No such 
structures have been observed in our simulations. Fur- 
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thermore, we observe some localized two-phase fingers 
that play a role during the instability of planar fronts at 
off-eutectic compositions, and that are similar to struc- 
tures seen in the experiments (compare, in particular, 
our Fig. |l2] to Fig. 6 of Ref. However, other exper- 

imentally observed patterns, such as "multiplet fingers" 
and two-phase dendrites are not reproduced by our simu- 
lations. It is possible that the existence of such patterns 
depends sensitively on the structure of the eutectic phase 
diagram, in particular on the asymmetry of the two solid 
phases and their surface energies that have been shown 
to influence the stability of binary lamellar eutectics , 
and on crystalline anisotropy. 

The present phase-field model could easily be modified 
to include some degree of asymmetry between phases as 
well as both solid-liquid and solid-solid anisotropy. In 
addition, the use of more general phase-field models with 
several order parameters j32|-[35[ , as well as the use of 
more efficient phase-field formulations [p9| and numeri- 
cal algorithms |l6 47 that greatly enhance the accessible 
length and time scales, could help to elucidate these ques- 
tions in the future. The exploration of the enormously 
vast parameter space of growth conditions and material 
properties that govern the formation of complex two- 
phase microstructures remains, however, a formidable 
numerical task. 
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